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IMAGE  UNDERSTANDING  AND  INFORMATION  EXTRACTION 
Research  Summary 

report  summarizes  our  research  progress  during  the  period  February 
1 *r..  30.  .976.  image  Understanding  and  information  Extraction. 

Our  research  objective  is  to  achieve  a better  understanding  of  image 
structure  and  to  use  this  kn^iedge  to  develop  technics  for  information  ex- 
traction from  imagery.  It  is  our  hope  that  the  results  of  this  research 
form  the  basis  for  the  deveiopment  of  technology  relevant  to  miiitary  appilca- 
ttons  of  machine  extraction  of  Information  fr»  aircraft  and  sateiiite  Imagery 
Our  research  projects  fa.,  into  five  heavliy  overiapping  areas:  image 

Segmentation,  image  Attributes,  image  Structure,  image  Recognition  Technigues. 
Preprocessing,  and  Applications. 

'mCE  SECMENTAT,0N  ' *"«  «•  approaches  to  iamge  segmentation: 

edge  detection,  and  region  growing.  In  edge  detection,  the  thrust  of  our 

research  is  to  use  syntactic  methods  to  heip  edge  detection  In  noisy  and  blur- 
red imagery.  As  a step  aiong  that  direction,  we  are  studying  characteristics 
of  digital  edges.  Aimost  a,,  past  works  on  digital  edges  and  curves  dealt 
with  the  Idea,  situation,  while  wear,  interested  mainiy  real-iife  digi- 
tized Images.  Our  experimental  results,  reported  by  Tang  and  Huane.  indicate 

that  the  properties  of  rea.-.ife  digital  straight  edges  are  quite  different 

from  those  of  Ideal  edn^c  in 

9es.  In  the  same  report,  a technique  of  recognizing 

real-line  digital  straight  edges  is  also  presented. 

A by-product  of  our  edge  detection  research  is  a technique  of  accurately 
estimating  edge  locations  which  holds  great  promise  in  mensuration  applica- 
tes. This  technique,  reported  by  Burnett  and  Huana.  is  based  on  the  Viterbi 

algorithm,  it  can  take  account  of  arbitrary  edge  profiles  and  fi,m  noise  and 
IS  computationally  slmpleo 
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In  region  growing,  we  are  studying  both  region  merging  techniques  and 

different  similarity  measures  among  regions. 

IMAGE  ATTRIBUTES  - We  are  doing  both  texture  and  shape  analysis.  A new 
class  of  texture  descriptors,  the  max-mln  descriptors,  has  been  developed. 

These  descriptors  are  very  simple  to  compute  and  perform  extremely  well  In 
various  classification  problems.  We  are  currently  exploiting  the  use  of  max- 
mln  descriptors  In  texture  boundary  detection  and  region  growing. 

In  shape  analysis,  we  have  done  extensive  study  on  Fourier  boundary 
descriptors  and  are  striving  to  Improve  their  performance.  We  are  also  devel- 
oping grammars  for  various  classes  of  objects  of  military  significance,  such 

as  airports,  and  tanks. 

IMAGE  STRUCTURE  - The  thrust  of  our  research  In  this  area  Is  to  use 
syntactic  methods  to  do  scene  analysis.  Fu  and.jj  report  on  the  us.  of  tree 
grammar  In  detecting  highways  and  rivers  In  LANDSAT  Imagery.  Tree  grammar  was 

also  used  In  helping  scene  segmentation  (Fu  and  Rena). 

IMAGE  RECOGNITION  TECHNIQUES  - We  pursue  two  topics:  the  use  b anc 

and  bound  techniques  In  solving  recognition  problems,  and  the  use  of  context 
In  statistical  classification.  In  Swain  and  report,  we  see  that  the  use 

of  context  Indeed  Increases  the  classification  accuracy.  However,  It  Is 
computationally  tedious  on  conventional  serial  computer.  W.  are  therefore 

looking  Into  the  use  of  the  llllac  A (via  the  ARPANET). 

PREPROCESSING  - The  aim  of  preprocessing  is  to  change  the  Image  to  a 
form  which  Is  more  convenient  for  Information  extraction.  Most  Images  are 
degraded  by  noise  and  blurring.  The  reduction  of  noise  and  the  sharpening  of 
the  Image  generally  facilitate  Information  extraction.  Several  noise  reduc- 
tion techniques  are  compared  by  Ygp  apd  puma.  It  was  found  that  the 
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synthetic  highs  technique  works  best  in  balancing  the  noise  level  and  the  edge 
sharpness 

APPL-\  , ONS  - We  are  working  on  several  applications,  two  of  which  are 
reported  here.  Wal lace  and  Wintz  present  some  recent  results  on  using  Fourier 
descriptors  for  airplane  classification.  M i tche 1 1 and  Chen  show  results  of 
reducing  cloud  and  haze  in  LANOSAT  imagery  using  three-dimensional  digital 
f i 1 ters. 
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J.  Burnett  and  T.  S.  Huang 

'■  i»  “ •' 

-» • 

, .-asurements  from  noisy  and  blurred  Images.  Here 

::;r: i.—— — - 

^A^»r\u»l.tl«.  of  a specific  Photographic  Imaging  systeo  slu«s 
the  algorithm  to  produce  asymptotically  efficient » and  unbiased  estimates  of 

object  sizes. 

, a Posteriori  Probability  «-T“"r«  Estimation 

S22*222 habl  1 1 tv  (MAP)  estimate  of  a sequence  i 

The  maximum  a-post«rlor  pro  deflned  ,5 

re  /_  belng  a degraded  and  noisy  version  o J 
given  a sequen  t p(jJi)  * is  a maximum.  To  calculate 

a sequence  ! - 0,.l2*--V  We  assume  the 

“1  a model  Is  needed  for  the  relationship  between  i _• 

- , /,  , I ) is  the  sequence  of 

observation  model  shown  In  Fig.  «•  1’  2-  ' ■ 

, , .,h  , the  nght  intensity  of  the  kth  samp,  h 

Ideal  light  intensities  with  k DOSSlble  values 

* m Each  i can  assume  one  of  G possible 
entering  the  Imaging  system.  k 

a a.  For  example  X might  represent  the  sequence  of  reflect 

al . aerial  photo  of  a bridge  across  a river. 

Intensities  from  a scan  line  o an  corresponding 

.a  b.  two  possible  Intensity  levels:  a,  correspon 

In  this  case  there  would  be  two  poss. 

. _n£j  g corresponding  to  the  g 

to  the  light  reflected  from  the  water  and  a2c 

— -rr. ::::: 

bridge).  The  state  at  position  k,  nk» 

, ..  , . ).  Since  each  I can  assume  only  a f.n. 

intensities  ( • k-v»  * * * * k k+v  J 
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number  of  values  each  „k  is  one  of  a finite  set  Is,.  . . s ].  Further  (t0 

within  boundary  conditions)  there  is  a one  to  one  correspondence  between  the 
state  sequence  r,  and  th»  Intensity  sequence  I. 

The  system  h(-)  represents  the  degradation  of  the  sequence  L the 
case  of  photographic  imagery  this  includes  blurring  due  to  scattering,  dif- 
fraction, camera  motion,  etc.  as  well  as  the  noniinear  relationship  between 

"9ht  3"d  fi’"  "»  oAiy  assumption  that  we  make  on  h is 

that  there  is  a one  to  one  correspondence  between  V = (V| ) (where 

yk  = h(n^))  and  n_„ 

- ‘S  3 SeqUa"Ce  °f  '"^Pendent  noise  samples.  We  do  not  rule  out 
dependence  of  the  noise  parameters  on  the  signal,  however,  Fbr  examp, e f„m 

grain  noise  is  approximately  norma!  with  . standard  deviation  proportional 
to  the  signal  level. 

The  Algorithm 

By  definition  the  MAP  sequence  estimate  I of  I is 


(i ) P(jJZ)  is  a maxi 


i mum 


bu,  since  there  is  a one  one  correspondence  between  I and^fl)  ,s 
equivalent  to 


(2)  P(njz)  ln=^  is  a maximum 
However, 

(3)  p fill Z.)  ■ P(Z,n)/pU)  - P rj_)  P(tl)/p(z) 

Due  to  the  Harkov  assumption  of  a and  the  independence  of  the  noise 

p(n>  -j!  P(nwlnk) 
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PCzffl)  - P(ztl) 


n h 

w p(zJYJ  * v 

k»l  * k k»l  k k 


Thus  we  want  to  maximize 

W » p<"w,ln,)  PUtMnJ) 


or  equivalent  minimize 


M „ 

£ - *»  p(vJV  ' £n  p^Jtih^t>)  “£z  r(V 

&y  assigning  a cost  or  length  of  F(hi^)  to  each  branch  of  a trellis  we  can 

A- 

see  that  the  MAP  estimate  ij  represents  the  lowest  cost  or  minimum  length  path 

through  the  trellis. 

£ 

Let  hi j be  a sequence  of  states  starting  at  some  initial  state  at  posi- 
tion one  and  ending  at  state  n at  position  £.  In  general  there  will  be 

SL 

several  possible  paths  (or  sequences  hi|)  through  the  trellis  that  pass 

through  state  hi  at  position  I.  Denote  this  set  of  paths  by  5.  Let  n*  e 5 

be  one  of  these  sequences  but  with  the  additional  restriction  that 
£ ~ 


F(hu)  m £ r(hi;)  < r(hi.)  for  any  other  sequence  hU£  H • (If  the  two  paths 
1 j-l  J 1 1 

have  equally  low  cost  any  reasonable  procedure  for  deciding  between  them 
will  do.) 

■*£ 

Thus  hi|  (called  the  survivor  sequence  or  survivor)  is  the  minimum  cost 
sequence  or  path  from  the  fixed  initial  state  to  the  state  ij  at  position  £. 


Now  if  the  minimum  cost  complete  path  from  position  I to  position  m passes 
through  state  hi  at  position  £ it  must  have  hij  as  its  initial  segment  (If  it 
did  not  we  would  replace  the  initial  segment  with  irji|  and  get  an  even  lower 
cost  path,  a contradiction). 


w.rHKUr '*■■*" 


Of  course  we  do  not.  know  that  the  minimum  cost  complete  path  passes 
through  state  q at  position  f ,,  However  we  do  know  that  it  must  pass  through 
one  of  p possible  states  at  position  5,.  Thus  at  any  position  l we  only  need 

A y A y 

store  at  most  p survivor  sequences  n j"  2nd  their  costs  f(q').  To  get  to 
position  5>1  we  need  only  extend  all  position  f.  survivors  by  one  unit,  com- 
pute the  costs  of  the  possible  extensions  from 


J r(r,®+l)  = r(n*)  * r(V|) 


and  for  each  of  the  p possible  states  q^j  at  position  to-1  select  the  lowest 
cost  path  ending  in  that  state  as  the  position  iid-l  survivor.  In  summary: 


Storage: 


"1 


(position  index) 

(p  such  1 point  survivor  sequences) 


F(h|')  (costs  of  each  of  the  p survivor  sequences) 


initialization:  too 


nQ  = qig  for  each  possible  initial  state 


F(qi l)  = -ton:  where  w is  the  a-priori  probability 


qi, 


•o 


that  the  initial  state  is  q_  (if  known) 

0 

if  the  a priori  probabilities  are  not  known  then  any  reasonable  initial 
cost  assignment  (such  as  F(q^)*0  for  all  possible  initial  states  q)  will  do. 
Recursion 

For  each  of  the  p possible  states  at  position  HL+l  compute 

r^rV  * r(nf)  + 

for  each  possible  q,r  find  r(n^+l)  = min  F(q; qj  store  F(qf+l)  and  the 
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corresponding  survivor  sequence.  At  position  M there  will  be  at  most  p 
survivor  sequences,  one  survivor  sequence  terminating  at  each  permissible 
position  M state.  Denote  by  n the  lowest  cost  of  these  sequences.  ^ Is  the 
MAP  sequence  estimate  of  rj. 

Example.  Suppose  it  Is  known  that  at  sample  point  number  one  the  local 

gray  level  is  zero,  that  at  position  eight  the  local  gray  level  Is  three  and 

that  someplace  between  these  two  points  a step  change  between  the  levels  zero 

and  three  occurred.  The  degrading  system  is  linear  and  shift  invariant  with 

impulse  response  h-l-h  -h  -i.  The  noise  is  white,  Gaussian  with  variance 
2 

a . The  observed  sequence  Z « (-.5,  +.25,  +.75,  .5,  2.8,  2.7,  3.3,  3.1). 

The  possible  states  are 

S,  “ S3  - (0,3,3) 

S2  " (0'°'3>  - (3,3,3) 

h(sj)  - 0 h(S3)  = 2 

h(s2)  - 1 h(S4)  - 3 

y assuming  that  all  permissible  changes  of  state  are  equally  likely 

"4n  P*nk+JV  terms  are  the  samc  and  can  be  ignored.  Thus  T(n  ) - (z0-h(n))2. 
The  trellis  for  this  example  is  shown  in  Fig.  2.  The  permissible  paths  are 
shown  by  dashed  lines.  The  numbers  are  the  costs  of  the  survivor  sequences  to 
each  state.  The  MAP  estimate  is  shown  by  a solid  line.  Its  total  cost  is 

A 

1.955  and  corresponds  to  1 - (0,0,0,0,3,3,3,3). 

The  minimum  cost  or  minimum  length  path  through  a trellis  problem  and 
various  solution  have  been  around  for  some  time  [2],  The  algorithm  presented 
above  (commonly  called  the  Viterbi  algorithm  or  VA)  was  presented  by  Viterbi 
[9]  as  a technique  for  decoding  convolutional  codes.  The  algorithm  has  since 
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been  used  by  Forney  [3]  as  a solution  to  the  intersymbol  interference  problem. 
The  algorithm  also  has  applications  in  text  recognition  [A]  since  it  can 
exploit  Markov  dependence  in  English  text  and  can  be  applied  to  scene  analysis 
problems  [5]  due  to  Markov  relationships  between  objects  in  scenes. 

The  VA  can  be  used  to  make  measurements  of  objects  in  digitized  images 
providing  the  object  whose  length  or  width  is  to  be  measured  is  distinguished 
from  its  background  by  abrupt  changes  in  reflected  light  intensity  at  its 
edges.  A MAP  sequence  estimate  of  a scan  line  across  the  object  can  be  cal- 
culated and  the  edge  locations  relative  to  the  start  of  estimated  line  can  be 
subtracted  to  produce  an  estimate  of  the  size. 

In  the  next  sections  formulas  will  be  derived  for  predicting  the  per- 
formance of  the  VA  at  locating  edges  and  estimating  widths. 

Error  Analysis  for  Step  Edges 

In  the  previous  section  we  presented  the  VA  as  a possible  technique  for  maki 
measurements  from  digitized  images.  The  measurement  technique  consists  of 
making  a MAP  estimate  of  a scan  line  across  the  object  of  interest  and  sub- 
tracting the  position  of  the  boundary  points  of  the  estimate  of  the  scan  line 
to  get  an  estimate  of  the  width.  In  this  section  we  examine  the  accuracy  of 
this  technique.  In  particular  we  will  calculate  the  probability  of  mislocat- 
ing  a step  edge  by  |n | (n = +i ,+2,. . .)  points.  It  is  assumed  that  a change  in 
brightness  from  level  a^  at  some  initial  position  to  level  af:  position  m 
has  been  detected.  The  observed  signal  is  blurred  and  noisy  so  that  the 
location  of  the  step  change  between  levels  a^  and  a^  is  uncertain. 

Define  an  error  event  E by  the  conditions  i.  =i.  and  i.  . . , = i,  , . , 

n k k k+|n|+l  k+jn|+l 

but  ij  t i.  for  k+1  <_  j £,k+|n|.  If  p is  the  number  of  states  then  there  must 

/v 

Ini  + p - 2 { j n | >_  1)  state  disagreements  (see  Fig.  3 for  p=4,  n=-2).  If 

A 

is  the  output  vector  corresponding  to  and  is  the  vector  corresponding 
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to  the  true  sequence  ^ then  ^will  be  decided  over  ^ If 

k+ | n | +t  a k+ | n I +t 

(7)  £ An  P(z.|y.)  > Z An  P(z.|y.) 

j-k-t  J J j-k-t  J J 


where  t » C (P-I )/2]  - largest  integer  <_  (P-1 )/2 

A 

If  the  noise  is  normally  distributed  with  zero  mean  and  variance  q then  (7) 
becomes 


(8) 


k+|n|+t  a ~ k+|n|+l 

T.  (z.-y.)2  < z 
j-k-t  J J j-k-t 


(zj'yj’ 


or  ||r-  r||2  < ||r-  Yjl 

where  ^ 

(Vt”",yk+t+|n|) 
(Vt W|n|> 


2 

Zk+  n +t 


) 


The  vectors  Y?  Y"can  be  viewed  as  three  points  In  |n|+2t+l  dimensional 
space  that  define  a plane  (see  Fig.  5). 

Define  ^Z,Y>  - z\  - Z z.y.  - Y^Z,  Now  the  condition  for  deciding 

j jj 


Y becomes 


[y=7|J 


(the  reason  for  the  normalization  by  | | Y— Y | | will  be  obvious  shortly) 
1 


-K- 

I Y— Y I 


- <z;r-z*  < <y;y?z*  - <r*,Y-z^] 


A A 


* — - <z;y*  < <r;£-z*  - <£y*] 


-Y 


|-s> 


."tf  ftfrtyrj’w,  <iii»>r]«W(*M'^-If,i/ii»K«(ir 


l^rtniwhll 


— A 

I Y-Y I 


f*I»  ?i'>  - 't'X>  - <rt  y">  < <y;  Y-rs  - <z;  y">  - <y;  y">; 


[<y;  r-z"> 


<y;  y-  t **>  < <y;  y^r>  - <y;  ^r>] 


A A A 


j_Y"  - Y"|  | | |Y"  - Y'| 


A 

Thus  an  error  will  occur  in  deciding  between  Y"and  T whenever  the  distance  from 
£ to  Y^*  a long  the  ~ - axis  is  less  than  the  distance  from  Z'to  Y'along  the 

nr- r ii  " 

y"-  r 

same  axis.  Since  Zf»  Y%  N_"the  projection  of  Z'-  Y'on  * — is  determined 


ir-  r i 

by  the  projection  of  the  noise  N"  on  the  axis.  This  projection  is  a linear 
combination  of  normal  random  variables  with  zero  mean  and  variance  a2  and  hence 
is  also  a normal  random  variable  with  zero  mean  and  (due  to  the  normalization 
by  1 1 Y”Y 1 1 ) variance  a . Therefore  the  probability  that  Y*  will  be  decided  over 
the  correct  path  Yj*  is  the  probability  that  a normal  random  variable  exceeds 

A A 

half  the  distance  between  Y“  and  Y'„  With  T = IlY"-  Vl|  then 

n 1 11 

(10)  Prob  (decide  Y"  over  Y“)  = /”  -—1—  e_n/2 ^ 

n / 2 to 
2 

- • 

A A 

Now  Y - Y = (|  - I)  * h (assuming  a linear  and  shift  invariant  degradation) 

A 

I - I - (0,0,0, A, A, A, 0,0) 


where  A = + (a^  - a^) 
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Let  U,(K)  be  the  unit  step  response  of  h.  Then  ||T'-  V'll2  . Ily  . yl  1 2 

- £ (v,-;,)2  . t2 

j J J n 

2 ^ A 9 

" " -Z1t(,J-'j>  * h] 

J«1  J J 

2 k+|n|+t 

"A  2 [U.(j-K)  - U.O-Jnl-K)]2 

j-k-t  1 1 

(11)  » A2  ^ 1 ru.(j)  - U.  (j-|n|)]2 
j—t  1 


For  any  value  of  n there  are  two  possible  edge  locations  estimates  that 
are  |n|  points  In  error  corresponding  to  location  estimates  that  are  n points 
before  and  n points  after  the  correct  location.  Thus  the  probability  of  an 
error  event  for  step  edges  with  known  levels  Is 


02)  P(En)  - 24(^2.) 
where  Tn  Is  given  by  (II). 

Thus  far  we  have  only  considered  the  case  of  constant  noise  variance. 

However,  for  some  types  of  noise  such  as  film  grain  the  variance  of  the  noise 
varies  with  the  signal. 


Consider  again  the  problem  of  trying  to  locate  a step  edge  between  levels 

a,  and  a2  but  now  assume  that  the  noise  variance  of  the  tth  sample  o2  depends 

on  y^.  Equation  (7)  becomes  2 

k+|n|+t  k+InU* 

(13)  _Z  ^ -in  (y^oc_)  - (z.-y.jV.  > Z - in  SR  a 


jck-t 


J J 


j“k-t 


- (YVjlW 


or  In  an  obvious  notation 


i¥'.'V''  iw-r1  «*■  * vj  i i*.  VJ/pf-A'.V  T ’V>?»  iT>  r*  fT'1  W.* A tiW.  Kw  * WBW*«n 


I151MWWW''  'WC iVl-'i**  <nWt-,7ll>)nrjl 
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(1<0  I'-  Y' 


0 2 < 1 1 L'-  l'\ |tf2.  + c* 

y yy 


where  C „ = 2!  ftn  -ft-  . 

yy  yj 

A 

The  exact  probability  of  deciding  Y over  Y can  be  calculated  from  (1A) 
by  a procedure  given  in  [6]  though  it  is  very  difficult  and  will  not  be  con- 

A 

sidered  here.  However,  one  observation  should  be  made,  if  Y is  the  path 

A 

through  the  trellis  that  mislocates  the  edge  by  +n  points  and  Y-n  the  path  that 
mi s locates  the  edge  by  -n  points  then  in  general  CA  j CA  and 


2 

I liL  “ 411%  t \\L'~  • Therefore  the  probability  that  Y is  chosen 

yn  y-n 

A 

over  Y_  will  not  be  the  same  as  the  probability  that  Y is  chosen  over  Y. 

— -n  — 

This  lack  of  equality  between  Pr(n=a)  and  Pr(n=-a)  causes  the  random  variable 
n to  have  a nonzero  mean  value  which  introduces  a bias  in  the  edge  location 

estimate.  Since  En  = EaPr(n«a)  and  since  Pr(n=a)  decreases  with  increasing 

a 

signal  to  noise  ratio  (SNR)  the  bias  will  decrease  with  increasing  SNR. 
Further  the  decision  boundaries  (and  hence  the  bias)  among  the  various 
possible  paths  are  determined  by  the  variances  0^  which  in  turn  are  deter- 
mined by  the  possible  levels  a?  and  the  degrading  function  h.  Thus  the 
probability  of  mislocating  the  edge  by  n points  (and  hence  the  bias)  is 
independent  of  the  sarnole  point  at  which  the  edge  occurred. 

If  the  SNR  is  high  enough  the  bias  can  probably  be  ignored.  If  the  SNR 
is  low  the  bias  can  be  calculated  by  the  procedure  mentioned  earlier  or  found 
experimentally  by  computer  simulation. 

Extension  to  Pulses 


yny  y-ny 


In  measuring  the  size  of  an  object  two  edges  must  be  located.  If  the 
object  size  is  sufficiently  large  or  the  signal  to  noise  ratio  is  high  enough 
then  with  high  probability  the  minimum  cost  path  and  the  correct  path  through 


I 
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the  trellis  will  coincide  somewhere  inside  the  object. 

With  this  assumption  the  events  Ej  » first  edge  mislocated  by  n1  points 

and  E2  13  second  edge  mislocated  by  n2  points  are  independent  [7].  Thus  if 

n is  the  number  of  sample  points  difference  between  the  true  width  and  the 
w 

estimated  width 

(12)  "w=n2-n, 

P(n  = a)  = E P(n,  = a)  P(n.  = a - 6) 
w g L 

where  P(n2  = a)  and  P (n 1 — a - B)  can  be  calculated  from  (10).  The  sum  is  over 

all  possible  values  of  3 that  n2  can  assume  such  that  n2  “ ni  ° a* 

If  the  noise  varies  with  the  signal  level  n may  not  have  zero  mean.  As 

w 

in  the  previous  section  the  bias  can  probably  be  Ignored  if  the  SNR  is  high 
enough  or  calculated  theoretically  or  found  experimentally  if  necessary. 

Unknown  Levels 

In  the  previous  section  we  presented  formulas  for  the  probability  of  mis- 
locating  an  edge  by  |n|  points.  This  analysis  assumed  that  the  poss j . e levels 
and  a2  that  of  the  process  were  known.  In  this  section  the  probability  of 
error  is  calculated  assuming  that  the  levels  are  not  known  a-prlorl. 

A reasonable  course  of  action  in  the  case  of  unknown  levels  Is  to  obtain 
"training"  samples  of  the  gray  levels  characterizing  the  object  of  interest. 

A A 

These  samples  can  be  processed  in  some  fashion  to  produce  estimates  and  a2 

A A 

of  levels  a j and  a2,  respectively.  Suppose  a^  = a^  + and  a2  “ a2  + e2 

2 

where  and  e2  are  normal  random  variables  with  zero  mean  and  variance  cr^ 

n A A 

and  ct2,  respectively.  (If  a1  and  a2  are  large  sample  maximum  likelihood 
estimates  of  a1  and  a2  then  the  above  model  is  accurate  [8].) 


m ■*  fB w«rw« 


j. n»i'.u . i‘ i-’,Mv«Vr.''!r‘ ''..»■  ?. 11 1 iVWT!  nrHir  lifter W 


•wRvwrm^  cu*« 
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Let  the  correct  path  through  the.  trellis  by  Y and  the  estimate  be  Y_.  In 
this  case  the  correct  path  Is  the  one  that  places  the  edge  in  the  correct 
position  even  though  the  levels  (or  height)  of  the  step  may  be  incorrect.  The 
"signal  space"  diagram  of  Fig.  5 has  been  repeated  in  Fig.  6 with  the  addition 
of  a new  point  0.  0,  represents  the  output  sequence  corresponding  to  the  true 
edge  location  and  step  levels.  In  general,  6 will  not  lie  in  the  plane 
defined  by  Z,  Y,  X-  However ,^as  before,  only  the  projections  of  X " L and 
v - 7 whprp  7 = 0 + N on  will  have  any  effect  on  the  choice  between 

7 - - - - TTy^yTT 

X and  X* 

Now  Y-Y  = Ah  * (0,0, . . u j , 1 , 1,»  • • " ^ - l^a| 

n 1 1 s 


and  Y-0  = 
of  Y-0  on 


h 


* (^,.e 


X”X 

Ty^yTT 


e,,e2,e2,...,r.2).  Thus  the  projection 
will  be  a linear  combination  of  and  e ^ 


T* 

i 

Since  a 1 inear 


combination  of  normal  r.v’s  is  again  normal  the  only  information  needed  to 
characterize  l'  is  its  variance  (it  will  have  zero  mean  since  e,  and  e2  were 

assumed  to  have  zero  mean)  o\  which  can  be  calculated  from  knowledge  of 
2 2 

a,,  ^2’  V and  h* 

Again,  by  arguments  similar  to  those  in  the  previous  sections  the 
probability  of  deciding  Y over  Y is  the  probability  that  the  noise  alone,  the 

Y-Y  axis  exceeds  ^ - T'or  equivalently  probability  that  the  sum  of  two  normal 

2,2  , Tn 

random  variables  with  variances  a and  exceeds  2 


(13)  = Q( — r — j-)  and 

2 (a  + a J 


P(E  ) = 2Q( j 0 

n 2(a  + a ) 

n t 
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This  result  can  be  extended  to  pulse  width  accuracy  as  before.  If 

nw  “ ni  " n2  ls  the  number  of  sample  points  of  error  In  the  width  the  estimate 
then 

Pr[n  =a]  = l Pr[n  =0]  Pr(n  =a-3) 

w 0 ; 2 

S A 

where  Pr(n.-B)  ■ Q( * — 5— ).  Equation  (13)  shows  that  If  the  estimates  a. 

' 2(o  «,)  1 

A o C 2 

and  a2  are  such  that  a < < a then  the  lack  of  a-prlorl  knowledge  about  the 
possible  levels  will  have  very  little  effect  on  the  probability  of  mlslocating 
an  edge. 

Simulation  Results 

A pulse  of  width  thrity  sample  points  was  generated  and  blurred  by  a 
linear  shift-invariant  system  with  a Gaussian  shaped  impulse  response  with  a 
standard  deviation  of  one  sample  point.  This  blurred  pulse  was  transformed  by 
^k  * ^ *066  0°9  1 ^ 1 • 5.)  + .2  where  1^  is  the  k— sample  of  the  blurred 

pulse.  This  transformation  simulates  the  d-lcg  E curve  of  film.  Noise  of 
standard  deviation  .^Yk'  J was  added.  The  noisy  blurred  signal  was  then  pro- 
cessed to  produce  an  estimate  of  the  pulse  width.  Several  different  SNR's 
were  used  with  one  hundred  width  estimates  obtained  at  each  SNR.  The  results 
are  shown  Iri  Figs.  6 and  7.  The  bias  In  the  width  estimate  can  be  seen  to 
decrease  with  Increasing  SNR  as  expected.  The  variance  of  the  estimate  can  be 
seen  to  decrease  with  increasing  SNR  and  appears  to  asymptotically  approach 
the  Cramer-Rao  lower  bound. 
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Figure  k Probability  of  error  calculation 
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Figure  5 Probability  of  error  calculation  for  unknown  levels 


Figure  6 Accuracy  of  the  VA 


vs.  SNR  (in  dB) 


DIGITAL  STRAIGHT  EDGES 
G.Y.  Tang  and  T.S.  Huang 

Studio'.  of  t ho  digitization  of  n %l right  line  have  been  made  by  Rosenfeld 
(I),  Freeman  f?].  Horse  [3],  Gan  far  [**],  Brons  [5].  A set  of  rules  has  been 
established  to  govern  the  digitization  of  an  ideal  straight  line  which  can  be 
described  by  a first  order  mathematical  equation.  But,  in  the  real  world, 
most  edges  which  appear  to  be  straight  to  our  eyes  do  not  fall  into  this  cate- 
gory. On  the  contrary,  they  suffered  from  noise,  degradation,  blurring,  etc. 

In  this  report,  we  are  going  to  show  some  experimental  results  on  real  straight 
edges  and  to  compare  them  with  the  ideal  case.  A simple  testing  algorithm  is 
then  developed  to  determine  if  a real  edge  is  a straight  one  or  not.  The 
application  of  the  testing  algorithm  can  be  found  in  various  areas  such  as 
syntactic  pattern  recognition,  character  recognition,  scene  analysis  and  ef- 
ficient contour  coding,  etc.  Our  ultimate  goal  for  this  study  is  to  use  sy- 
tactic  method  to  aid  us  to  detect  edges  in  a noisy  environment. 

In  the  following,  we  are  to  use  the  word  'real1  to  refer  to  whatever  is 
on  pictures  obtained  by  practical  imaging  devices.  The  word  ‘ideal1  refers  to 
the  ideal  case. 

I . PROPERTIES  OF  AN  IDEAL  STRAIGHT  LINE: 

An  intensive  study  of  the  properties  that  a digitization  of  a straight 
line  should  have  has  been  reported  previously.  A brief  summary  is  given  here. 

(A)  The  Chord  property,  proposed  by  Rosenfeld,  is  a necessary  and  suffi- 
cient condition  for  a digital  arc  to  be  straight. 

(B)  The  digitization  of  any  real  curve  can  be  expressed  in  terms  of  a 
chain  code  [2],  The  chain  code  of  a straight  line  should  obey  the  following 
rules. 

1)  There  are  at  most  two  different  code  elements  differing  by  1,  modulo  8. 
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2)  One  of  these  two  codes  occurs  singly. 

3)  Successive  occurrences  of  the  single  element  are  as  uniformly  spaced 
as  possible. 

t)  For  an  ideal  straight  line  with  rational  slope,  the  corresponding  chain 
code  is  a repetition  of  a certain  period. 

D)  We  can  always  find  a rational  slope  to  approximate  a real  slope  so  that 
their  chain  codes  are  as  close  to  each  other  as  we  want. 

E)  An  ad  hoc  testing  algorithm  has  been  proposed  [3]  to  see  If  a given 
chain  code  may  come  from  a straight  line. 

II.  EXPERIMENT 

The  nice  properties  listed  in  the  foregoing  section  are  based  on  the 
assumption  that  the  lines  are  ideally  straight,  the  digitization  Is  noise-free 
and  the  digitization  scheme  Is  well  defined  [I],  [2].  |n  the  real  world,  where 
a model  with  such  a high  order  of  idealism  can  hardly  find  many  applications, 
minor  randomness  corrupts  the  structures  severely.  The  Intention  of  our  ex- 
periment Is  to  observe  the  code  structures  of  real  edges  and  to  compare  them 

to  the  Ideal  case.  Also  we  outline  a new  testing  algorithm  from  our  experi- 
mental results. 

A set  of  12  pictures  has  been  taken  by  a standard  Nikon  F-2  camera  with 
35  mm  film.  A flying-spot  scanner  Is  used  to  discretize  the  picture  Into  a 
square  matrix  of  size  256  by  256.  The  grey  levels  range  from  0 to  255.  A 
simple  contour  follower  and  a chain  code  encoder  Is  used  to  obtain  the  chain 

on  each  plcture‘  Only  straight  edges  are  on  these  12  pictures.  There 
are  totally  30  straight  edges. 

A contour  follower  is  designed  so  that  the  output  of  the  contour  follower 
is  the  chain  codes  of  the  edge  which  It  follows.  The  Input  to  the  contour 
follower  Is  a window  of  size  3x3.  We  assume  that  the  center  of  the  window  Is 
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on  the  contour  and  that  we  know  which  of  its  eight  neighbors  is  also  on  the 
contour.  The  task  of  the  contour  follower  is  to  find  a point  which  is  supposed 
to  be  the  next  contour  point  among  the  eight  neighbors  around  the  center.  Then 
the  window  moves  to  and  centers  at  that  point.  The  direction  of  move  is  en- 
coded by  Freeman's  method.  Repeating  the  same  procedure,  we  can  obtain  the 
whole  contour.  The  way  to  locate  the  next  point  from  the  eight  neighbors  of 
the  center  is  simply  looking  at  the  eight  neighbors  in  counter-clockwise  sense 
starting  with  the  neighbor  already  known  on  the  contour.  A thresholding  tech- 
nique is  used  to  discriminate  the  two  regions  defining  the  contour. 

Table  1 shows  which  ideal  properties  are  violated  by  chains  obtained  from 
real  pictures. 

The  ideal  edge  is  not  what  one  finds  in  the  images  produced  by  real-life 
imaging  devices.  There  are  several  factors  that  degrade  the  edges  that  are 
actually  found.  Two  predominant  ones  are  blurring,  or  defocusing,  and  ir- 
regularities of  the  surface  structure  of  the  object.  Besides,  for  the  case  of 
straight  edges,  a slight  concavity  or  convexity  can  be  considered  as  another 
kino  of  disturbance.  The  net  effect  of  these  disturbances  of  Freeman  chain 
code  can  be  summarized  as: 

1)  A third  code  element  is  introduced  if  there  is  a missing  or  spurious 
lattice  point.  This  third  code  element  will  occur  together  with  the  code 
element  which  occurs  singly.  Figure  1 illustrates  this. 

2)  A third  run  length  is  introduced  if  the  missing  or  spurious  lattice 
point  occurred  at  the  beginning  or  at  the  end  of  each  run.  Fig.  2 illustrates 
this. 

3)  Other  types  of  errors  occur  which  cannot  be  characterized  easily. 

Fig.  3 shows  us  the  test  pictures. 
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III.  STRAIGHT  LINE  RECOGNITION 

To  recognize  a real  chain  code  as  a straight  line  is  not  straightforward 
since  we  have  no  rigorous  definition  of  what  is  a real  straight  edge.  Here  we 
propose  to  employ  a heuristic  method  to  soive  the  problem  of  the  recognition 
of  digital  straight  lines. 

The  basic  strategy  of  our  approach  is  rather  simple.  We  are  attempting  to 
enclose  a digital  straight  line  by  a straight  band.  If  such  a band  with 
reasonable  bandwidth  could  be  found  for  a given  chain  code,  then  we  claim  that 
that  chain  code  corresponds  to  a straight  line.  Mathematically,  a straight 
band  is  defined  by  an  equality 

|Y  - mX|<_k 

where  m is  the  slope  of  the  straight  band  and  k is  related  to  the  bandwidth  of 
the  straight  band.  Two  parameters,  m and  k,  are  therefore  to  be  determined 
from  a given  chain  code.  We  proceed  as  follows: 

Step  1 « Check  if  there  are  only  two  code  elements;  one  occurs  singly. 

Yes;  go  to  Step  3 
No;  go  to  Step  2. 

Step  2 * Correct  all  possible  first  kind  of  errors  mentioned  in  the  pre- 
vious section.  Then  check  if  there  are  only  two  code  elements  left;  one 
occurs  singly.  If  it  is  yes,  then  go  to  Step  3;  otherwise  do  the  following; 

Calculate  the  percentage  of  the  3rd,  4th  ...  etc.  code  elements  to  the 
total  length  of  the  chain  code  and  the  percentage  of  the  second  code  elements 
which  should  occur  singly  in  ideal  case  but  occurs  in  runs  in  this  case  to  the 
total  number  of  the  second  code  elements.  Then  compare  these  two  percentages 
to  two  preset  thresholds.  Once  they  are  smaller  than  the  thresholds,  go  to 
Step  3;  otherwise  the  chain  code  is  rejected  as  a straight  line. 


% 

| 

■4 
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Step  3 ~ m i s the  ratio  of  the  element  which  occurs  singly  to  the  total 
number  of  1st  and  2nd  elements. 

k is  a subjective  quantity.  It  relates  to  the  noise  level  of  the  picture, 
the  sharpness  of  the  edge  and  the  resolution  of  the  picture. 

It  seems  that  so  far  we  have  taken  few  advantages  of  the  nice  structural 
regularities  existed  among  ideal  straight  lines  to  aid  us  in  developing  test- 
ing algorithm  for  the  non-ideal  case.  A rather  straightforward  heuristic 
method  has  been  devised  to  further  reduce  the  number  of  testing  points.  It 
is  noted  that  the  straight  band  is  a geometrically  convex  set.  So,  only  the 
two  ends  of  each  run  on  the  chain  codes  need  to  be  tested.  Furthermore,  if  we 
apply  a simple  operation  which  is  to  replace  each  run  by  its  run  length  and  to 
delete  the  code  elements  occurred  singly,  then  the  output  of  the  operation 
should  obey  the  rule:  only  two  code  elements;  one  occurs  singly,  for  the 
ideal  case.  For  the  real  chain  code,  we  may  apply  the  same  operation  itera- 
tively to  it  until  the  output  fails  to  follow  the  rule.  We  then  use  the 
final  one  as  a clue  to  break  the  initial  input  chain  into  several  pieces  of 

i 

line  segments.  Only  the  ends  of  each  line  segment  are  to  be  tested.  An 

example  illustrating  this  follows: 

A chain  code  from  picture  No.  7 in  our  experiment  is 

A 

010101010101010101010  i 00101010101010101 
A A 

0 10  10  10  10  10  0 10  10  10  10  10  10  10  10  10  0 1 0 10  10  1 0 

, A A 

101010101010100101010101010101010010101 

A 

01010101010101001010101010  1. 


30 


The  first  output  of  our  operation  is: 

A A A 

M M M » i I l 1 2 M I 11  1 ! I i I 1 1 2 11  1 I I I 1 1 2 

A A A 

I I II  i * M 1 1 2 1 11  11  1 1 1 2 1 1 1 11  I 1 1 1 2 1 11  1 1. 

The  next  output  is: 

^ *2  8 10  8 9 5 . It  fails  to  follow  the  rule.  So 

we  trace  back  and  locate  the  potential  testing  points  (wlthA  ) on  each  gen- 
eration. The  coordinates  of  the  testing  points  with  respect  to  the  initium 
are  then  (22,11),  (1*9,2/*),  (68,33),  (90,l*i») , (110,53),  (131,63). 

Table  2 shows  the  result  of  applying  the  foregoing  method  to  test  17 
chain  codes  obtained  in  our  experiment.  The  "MIN  THRESHOLD"  means  the  small- 
est k value  we  have  to  choose  in  order  to  report  the  edges  straight.  Also  we 
applied  Morses  algorithm  to  test  these  real  chain  code.  It  does  not  work  as  well. 
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Line  No,  I 

2-1  X 

2-2  X 

2-3  X 

2- 4  X 

3- 1  V 

3-2  V 

3-3  V 

3- 4  V 

4- 1  V 

4-2  V 

4-3  V 

4-4  V 

4-5  V 

4-6  V 

4-7  V 

4-8  V 

5 V 

6 X 

7 V 

8 V 

9-1  V 

9-2  V 

10-1  V 

10- 2  V 

11- 1  V 

11-2  V 

12  V 

1-1  X 

1-2  V 

1-3  X 


2 3 4 5 6 

x X X X X 

x X X X X 

x X X X X 

x X X X X 

V V V V x 

V X X X X 

v X X X X 

V X X X X 

V X X X X 

V V V V X 

V V V X X 

V V V V X 

V X X X X 

V X X X X 

V V v v V 

V V V V v 

V X X X X 

X X X X X 

V V V V X 

v X X X X 

V V V X X 

V X X X X 

X X X X X 

V X X X X 

v X X X X 

V X X X X 

V X X X X 

X X X X X 

V V V V v 

X X X X X 


1-  Only  two  code  elements 

2.  One  of  the  two  code  elements  occurs  singly 

3.  Only  two  run  lengths 

4.  Two  run  lengths  are  consecutive  integers 

5.  One  of  the  two  run  lengths  occurs  singly. 

6.  The  runs  of  the  run  lengths  can  be  of  at  most  two 

run  lengths:  one  of  which  occurs  singly 


Table  1,  Comparison  of  real  cases  and  ideal  case 
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Line  No. 

Length 

Min.  Threshold 

HORSE 

3-1 

147 

.795 

NO 

3-2 

49 

.755 

NO 

3-3 

54 

CO 

CO 

oo 

• 

NO 

3-4 

58 

.724 

NO 

4-1 

82 

1.329 

NO 

4-2 

47 

.404 

YES 

4-3 

42 

1.524 

NO 

4-4 

127 

1.055 

NO 

5 

238 

1.345 

NO 

6 

237 

1.456 

NO 

7 

130 

.435 

YES 

8 

209 

1.004 

NO 

9-1 

85 

M7 

NO 

9-2 

245 

.942 

NO 

10-1 

147 

1.258 

NO 

10-2 

40 

.900 

NO 

11-2 

181 

1.127 

NO 

Tab,e  2 inszii  rc-^rse4 
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perfect  case 


X - X - X - X - X - X - X - X - X 


X - X - X - X - X - X - X - X - X 

00000000  1 oooooooo 


a missing  point: 


x-x-x-x-x-x-x-x 


X - X - X X - X - X - X 

X 

0071000100000000 

two  missing  point: 


X - X - X 


x-x-x-x-x-x-x-x 


X - X - X 


X - X 

00701001000000 

a spurious  point: 

X x-x-x-x-x 

x - X - X X X - X - X - X 

0017000100000000 

two  spurious  point: 


- X - X - X 


X - X 


X - X - X X X X - X - X 


X-X-X-X-X-X-X-X 


00  10700100000000 


fig.  lu  Missing  or  spurious  points  occurred  in  runs 


3% 


perfect  cme: 

W ^ M 

^ ^ n' 

U«  Mr  Mr 

A*  A* 

W _ W _ W v W 

A A A A 

GDG)©  T1  Q)©©  H ©©(D) 

M MT  Mr  u«r 

A>  A'  A' 

a spurfioiuis  poiimtt: 

Mr  (W  \W  (W  M 

A«  — A—  A.-Xv-A 

W VK  \w  W 

A • A ^ A * A 

®©T1  ©®©©T1©©© 
a imiissiinDg;  poiiimts 

M »W  Mr  M> 

rv  fv 

M1  M \W 

w /v  A 

W \W  Ut*  M VW 

^ ^ /V  ^ ^ 

©©©©n  ©©d©©© 

Fiig.  2 Mtiissiirog  or  sporiiou®  poirots  ©c cored  at  tttfre 
oiMtis  of  a rum 
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No.  2 


No.  4 


Figure  3 Examples  of  test  pictures 
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No.  7 


Figure  3 Examples  of  test  pictures 


K.  S.  Fu  and  R.  Y.  Li 

The  essential  problem  in  our  research  on  the  syntactic  pattern  recognition 
of  highways  and  rivers  is  really  to  find  a grammar  that  will  describe  well 
these  classes  of  interest.  If  the  physical  shape  of  the  class  under  consider- 


ation is  completely  known  and  fixed,  like  printed  English  character;  we  can 
iiwnediate’y  write  down  the  syntactic  rules  to  describe  its  structure.  Since 

i 

this  is  not  quite  true  in  our  case,  the  construction  of  grammatical  rules  has 
to  be  based  on  informations  obtained  from  a set  of  sample  patterns  known  to 
come  from  that  class.  Hopefully,  this  set  of  inferred  rules  should  be  able  to 

I 

describe  and  predict  other  sample  patterns  which  are  of  the  similar  nature  as 
the  original  training  samples  and  presumably  in  the  same  class.  A basic  ap- 
proach of  grammatical  inferrence  problem  is  to  construct  a grammar  by  identify- 
ing the  syntactic  structures  of  the  known  string  and  any  possible  recursiveness 


that  might  happen.  There  are  three  steps: 

(1)  Try  to  discover  the  syntactic  structure  of  the  given  string  by  looking  for 
repetition  and  dependent  relationships. 

(2)  Decide  what  sublanguages  make  up  the  language  and  generate  non-terminals 


for  e.i.-h  sublanguages. 

(3)  Combine  equivalent  nonterminals  which  have  almost  the  same  sublanguage  and 
determine  the  appropriate  relationships  among  sublanguages. 

One  practical  method  to  learn  the  syntactic  structures  of  the  given  pictures 
is  to  use  a semantic  teacher  to  learn  the  meaningful  nonterminals  one  level  at 
a time  [1].  To  start  the  inferrence  process,  we  first  find  the  types  of  term- 
inals or  primitives  that  will  fit  the  subparts  of  the  picture  pattern  for  a 
given  window  size.  After  this  initial  extraction  process,  we  have  to  decide 
the  most  probable  combinations  of  primitives  which  occur  as  neighbors  of  each 

1 

} 

j 
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other  In  the  set  of  observed  training  samples.  These  combinations  are  then 
applied  to  the  training  data  set  to  test  their  recognition  effectiveness.  When 
the  results  appear  to  be  satisfactory  after  some  additions  and  deletions  of  the 
combination  rules,  we  can  choose  this  set  of  rules  to  represent  the  training 
samples.  The  appropriate  grammar  can  then  be  formulated  by  using  these  rules. 
In  our  present  case,  we  choose  the  tree  grammar  because  of  Its  easiness  to 
describe  these  rules.  A tree  recognition  program  based  on  this  tree  grammar 
can  then  be  used  to  recognize  the  training  data  set  of  Lafayette  and  a test 
data  set,  that  of  Grand  Rapids,  Michigan.  Figures  (1),  (2),  and  (3)  contain 
the  results  from  Lafayette  experiment.  Other  preliminary  results  from  Grand 
Rapids  can  be  found  in  reference  [3]* 
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Figure  (l)  Statistically-classified  results  of  Lafayette 
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Figure  (2)  Skeltons  from  Lafayette  samples 
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Figure  (3)  Syntactically-recognized  results 
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SYNTACTIC  SCENE  SEGMENTATION 
K.  S.  Fu  and  J.  Keng 

A syntactic  approach  to  scene  segmentation  has  been  investigates  which 
involves  two  levels  of  processing.  The  first  level,  referred  to  as  the 
transformation  process,  consists  of  five  steps  referred  to  as  (1)  threshold 
finding,  (2)  horizontal  processing,  (3)  vertical  processing,  (4)  logic 
integrating,  and  (5)  line  smoothing.  The  second  level,  which  is  the  actual 
syntactic  analysis,  requires  inference  of  a tree  grammar  to  describe  the 
boundaries  of  homogeneous  regions.  The  tree  grammar  is  then  implemented  in  a 
parser  which  traces  the  region  boundaries. 

The  approach  has  been  Implemented  and  initial  experiments  on  mu lti- 
spectral  remote  sensing  imagery  are  being  conducted.  Further  detail  and 
experimental  results  will  follow. 

Working  on  the  problem  of  picture  segmentation  through  a syntactic 
approach,  we  feel  that  the  evaluation  of  the  earth  resources  is  very  useful. 

So,  the  multi  spectral  remotely  sensed  picture  is  chosen  as  data. 

There  are  two  levels  of  processing  for  the  syntactic  picture  segmenta- 
tion, first,  the  transformation  process  and,  second,  the  tree  grammar  analysis. 
The  first  level  consists  of  five  sub-processes,  threshold  finding,  horizontal 
processing,  vertical  processing,  logic  integrating,  line  smoothing.  The 
process  of  tree  grammar  analysis  utilizes  the  correspond  .g  parser  from  the 
inferred  tree  grammars  to  process  the  transformed  picture.  Then  a picture  is 

segmented. 


• flr”  '—I  - t ransformation  processing 

(I)  Threshold  finding 

From  a digitized  picture  • 

’ draining  regions  are  located  fr 

homogeneous  parts  of  eh*  • * d from  eVery 

- - - r - — - - 

they  ere  close  t ' "" *"  If 

selected  es  a thresZT  ^ ^ ^ ““ 

(2)  horizontal  processing 

The  model  of  the  multisectral 

n-dimensiona,  space  E*  , h *•  Euclidean 

x x " ,S  3 $Pace  havi"9  n coordinates 

I’  2’“”V  The  number  n here  represents  the  h 

he  thosen.  A point  of  th  °f  Cha""el5  t0 

■—  - by  definition  an  ordered  „-tuple 

Eucnde”’’d-"  ' ^ d'SCanCe  bCtM!en  POi"tS  is  *«"*  as 

ean  dl stance.  The  operation  procedure  is  to  co 

levels  of  (l,|)  and  0,2).  ,fthetl,  . *"* 

then  set  zeroes  to  (|  ,)  and  (,  Sma,,er  ^ thresho'd, 

•f  the  distance  is  ^ «•»  «-  0.3,  are  chared. 

greater  than  threshoid.  A one  is  put  to  (,  „ , d 
the  same  operations  start  from  (,  M .. 

he  sa,,  pattern. 

(3)  Vertical  processing 

The  operations  are  the  . . 

the  same  as  horizontal  processing  except 

it  goes  verfiratiw  • y C*cept 

3 vc‘ti  ca  l l y instead  • 

1 r,b ceaa  of  horizontal  ly. 

(^)  Logic  integrating 

The  logic  variable  of  horizontal  procession  • 

Tor  the  vert.V.l  'S  named  as  H and  V 

vertical  processing^  The  iogic  integral 

the  integration  th  9 process  achieves 

9 at,°n  thr°U*h  B°°1ean  algebra  V+H. 
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(5)  Line  smoothing 

The  line  smoothing  algorithm  connects  the  discontinuity  of  the 
lines.  As  for  this  part,  lots  of  line  smoothing  algorithms  can  be 
devised. 

b.  Second  level  - tree  grammar  analysis 

A tree  grammar  is  inferred  to  describe  the  boundaries  of  the 
homogeneous  regions.  The  tree  grammar  is  used  to  trace  the  boundaries 
and  reject  the  unnecessary  boundary  parts.  Finally  a segmented  picture 
is  received  from  the  (two  level)  syntactic  method. 

The  scheme  of  syntactic  picture  segmentation  has  been  implemented  and 
the  experiments  have  been  conducted  on  the  mul ti spectral  remotely  sensed 
pictures  of  an  area  in  the  State  of  Indiana.  The  pictures  were  taken  on 
August  13,  1971  and  stored  in  the  computer  IBM  360/67  of  the  Laboratory  for 
Applications  of  demote  Sensing  in  West  Lafayette,  Indiana. 

The  result  of  a picture  96x96  area  is  shown  in  Fig.  1.  For  the  purpose 
of  comparing  processing  time  and  accuracy,  a result  of  statistical  segmenta- 
tion by  clustering  is  also  provided  in  Fig.  2.  The  computer  processing  time 
on  IBM  360/67  for  the  same  area  of  the  picture,  the  syntactic  picture  seg- 
mentation takes  about  36  seconds  and  the  clustering  technique  requires  180 

seconds. 

The  area  of  column  (105-132)  and  row  (444-450)  in  Fig.  1,  shows  lots  of 
boundary  points.  From  a survey  of  the  ground  truth  (Fig.  3).  it  points  out 
the  reason.  Because  this  area  is  pasture  and  bare  soil  compound  area.  So 
lots  of  boundary  between  these  two  objects  of  the  compound  area  are  located. 
A classified  result  of  the  same  area  is  also  provided  in  Fig.  4.  In  the 
region  column  (90-100),  row  (480-488)  in  the  syntactic  segmentation  result 


I 

V' ' 


*MK?i  W«v  ws-o* 


^5 


'V*  « ’’;•*, ^>*3131^3 


shows  a segment  which  corresponds  to  the  same  location  in  the  classified 

I 

result.  But  the  clustering  result  can  not  segment  it  well. 
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THE  SYNTACTIC  PICTURE  SEGMENTATION  RESULT 
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Figure  ] The  syntactic  picture  segmentation  result 
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Copy  available  to  DDC  does  not 
permit  iully  legible  reproduction 
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Figure  2 The  statistical  picture  segmentation  Clustering) 
on  same  area  as  Fig.  1 
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STATISTICAL  DEPENDENCY  MODELS  OF  CONTEXT 
P.  H.  Swain  and  T.  S.  Yu 

In  the  previous  quarterly  report,  we  described  a model  for  incorporating 
context  in  the  image  analysis  process  through  compound  decision  theory. 

Briefly,  under  a fairly  stringent  set  of  assumptions  a Bayesian  strategy  is 
employed  which  classifies  a point  into  one  of  a candidate  set  of  classes  based 
on  the  multispectral  data  from  the  point  itself  and  the  data  from  either  the 
neighboring  four  or  neighboring  eight  points.  Figure  1 shows  the  results  of 
applying  this  approach  for  classifying  a small  set  of  LANDSAT  multispectral 
scanner  data.  A block  of  imagery  128x128  pixels  was  classified  using  the 
"simple1"  rule  (no  context),  with  the  ^-neighbor  rule,  and  with  the  8-neighbor 
rule.  Samples  of  900  pixels  (30x30)  were  selected  exhibiting  a range  of 
pointwise  accuries  and  the  corresponding  accuracies  obtained  using  the  context- 
incorporating  rules  were  tabulated  together  with  the  corresponding  classifica- 
tion times.  The  results  of  the  experiment  show,  as  expected,  the  classifier 
using  context  is  consistently  better  than  the  classifier  which  makes  each 
decision  based  on  data  from  the  individual  points.  Furthermore,  the  "8- 
neighbor  classifier"  is  consistently  better  than  the  "^-neighbor  classifier". 

However,  the  price  paid  in  terms  of  computation  time  is  substantial.  To 
justify  general  use  of  this  approach  we  shall  have  to  demonstrate  (a)  its 
performance  potential  over  a sufficiently  broad  range  of  analysis  problems, 
and  (b)  a means  of  implementation  (special  purpose  hardware,  software,  or  a 
combination  thereof)  which  is  efficient  enough  to  provide  results  on  a cost- 


effective  basis. 


k neighbor  rule 


Fig.  1 Results  of  an  experiment  in  utilizing 
statistical  dependencies 
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IMAGE  NOISE  REDUCTION 
M.  Y.  Yoo  and  T.  S.  Huang 

I.  INTRODUCTION 

There  have  been  many  attempts  to  enhance  images  degraded  by  detector  noise 
and  imperfection  of  the  imaging  system  [1-4] . In  most  cases  we  encounter  two 
contradicting  requirements:  reducing  the  noise  as  much  as  possible,  and  retain- 
ing edge  sharpness.  The  only  reasonable  answer  is  a compromise  between  the 
two.  We  will  approach  this  problem  using  the  synthetic  highs  technique  where 
we  can  treat  the  smoothly  varying  part  and  the  sharply  changing  boundaries 
separately  and  we  may  enjoy  some  freedom  in  compromising  the  two  situations 
[5].  The  performances  of  the  system  will  be  compared  with  several  available 
heuristic  approaches  [6] t [7]. 

II.  IMAGE  ENHANCEMENT  SYSTEMS 

2. 1 Synthetic  Highs  Technique 

This  technique  was  proposed  by  Schreiber  [8]  and  was  used  in  two  dimen- 
sional contour  coding  for  data  compression  by  Graham  [9].  The  basic  idea  of 
this  technique  is  to  decompose  the  image  into  two  major  elements  (slowly 
varying  "lows"  and  synthetic  highs"  which  are  mostly  boundaries  and  textured 
parts)  such  that  the  recombination  of  the  two  elements  results  in  the  original 
again. 

The  technique  was  described  in  detail  in  previous  reports  and  in  Graham 
[9],  and  we  are  not  going  to  repeat  that  here  again.  To  use  the  technique  for 
noise  reduction,  we  need  a noise  reduction  filter  between  the  edge  detector 
and  the  reconstruction  filter.  So  we  need  an  acceptable  boundary  detector  for 
noisy  images.  Tang's  [10]  or  the  following  simple  edge  detector  may  be 
used : 


p(x,y) 


Figure  1 Synthetic  highs  technique 

First  we  calculate  the  grey  level  differences  in  four  different  direc- 
tions; horizontal,  vertical,  diagonal  1 (1*5°),  and  diagonal  2 (135°)  as 


Horz  = iG^rAyj-G  jx^y^Yjl  . lG(x-Ax,y)  - G (x+Ax.v) 

GU.y-AyJ  + G ix,y+Ay)  * e 1 Glx-Ax,y)  + 6(x+Ax,y) 


diag  1 = iiU-Ax>y+Ay)-G(x+Ax,y-Ay) 
G (x-Ax,y+Ay)  + G(x+Ax,y-Ay) 


diag  2 


G(x-Ax.y-Ay)  - G(x+Ax,y+Ay 
G(x-Ax.y-Ay)  + G(x+Ax,y+Ay 


Whenever  one  of  the  following  conditions  holds  we  decide  G(x,y)  is  on  the 
boundary. 


i)  Horz  6j  and  Vert  < 

ii)  Horz  < ©2  and  Vert  >_  0^ 

iii)  diag  1 >_  0^  and  diag  2 < 6^ 

iv)  diag  2 < 6^  and  diag  2 > 0^ 

where  0^  ©2  ( 0 <_  0Jf  0£  < 1)  are  preassigned  threshold  values.  Typical 
values  of  6| , ©2  are  0.2,  0.1. 
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2.2  Median  Filter 

We  use  an  nxn  (n  is  odd)  window  to  scan  the  noisy  image  and  replace  the 
grey  level  of  the  center  by  the  median  grey  level  within  the  window.  Pratt 
[6]  used  one-diniensionai  window  to  remove  impulse  noise.  Resolution  of  the 
filtered  image  is  highly  dependent  upon  the  size  of  the  window  and  appropriate 
size  should  be  chosen  to  retain  reasonable  boundaries.  3x3  and  5x5  windows 
were  used  for  our  experiment. 

2.3  Variable  Width  Filter 

First  we  take  the  gradient  of  the  noisy  image  and  divide  the  absolute 
gradient  level  into  four  different  intervals.  We  assign  the  smoothing  filter 
of  an  appropriate  size  to  each  interval  such  that  the  lower  the  absolute 
gradient  level  is,  the  larger  the  duration  of  the  filters  impulse  response  is. 
One  dimensional  Gaussian  and  median  fi Iters  can  be  used  in  the  horizontal  and 
the  vertical  directions  sequentially. 

The  basic  idea  of  this  approach  is  that  we  retain  more  boundaries  where 
the  absolute  gradient  level  is  high  by  using  narrow  smoothing  filters  and 
heavily  smooth  out  slowly  varying  part  by  wide  smoothing  filters.  Uniform  or 
non-uniform  subdivision  of  gradient  levels  may  be  used  depending  upon  the 
distribution  of  gradient  levels.  The  sizes  of  the  filters  used  areO,  3,  5,  7. 

2.4  Noise  Cheating  Technique 

This  technique  is  a combination  of  two  averagings  with  different  window 
sizes.  We  average  the  noisy  image  using  an  mxm  window  and  average  the  noisy 
picture  again  by  an  n x n (n  > m)  window.  In  the  original  paper  [73  the 
authors  quantize  grey  levels  of  "severely"  averaged  images  using  a quantum 
step  that  is  at  least  four  times  the  standard  deviation  of  the  averaged 
picture.  But  if  the  standard  deviation  of  the  averaged  picture  is  large 
enough,  the  quantizing  process  wipes  out  everything  and  this  really  happened 
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In  0ur  case,  so  we  may  have  to  .kip  the  guantiting  step  depending  upon  the 
sjze  of  the  standard  deviation.  After  averaging  the  noisy  Image  -■«-  *» 
different  size5  of  windows  we  rep, ace  each  grey  ieve,  in  the  -ii.hk.y  averaged 
image  hy  the  ciosest  grey  ieve,  aTOng  the  corresponding  eight  surrounding 

image  points  in  the  "severely"  averaged  image. 

i . • i-he  "1  iqhtly"  averaged 
The  idea  is  that  we  retain  the  resolution  of  the  light 

picture,  while  we  en,oy  the  reduction  of  noise  level  of  the  "severely"  averaged 

iraa9e.  The  combining  sch^e  is  a hind  of  discrete  maximum  iiheiihood  approach. 

2.  n = 3 were  used  for  the  experiment. 

m fvpfr 1 MENTAL  RfSiilTS  AND  CONCLUSION 

” white  Gaussian  noise  was  added  to  generate  a ,0  dB  (variance  of  signa,/ 

variance  of  noise,  noisy  picture  of  site  256x256.  The  bandwidth  of  the  low 

. n ii/i  The  original  noise~ 
pass  Gaussian  filter  in  synthetic  highs  system  ,s  0.1,6. 

• Fin  2 and  the  10  dB  noisy  picture  an- 
less  picture  is  shown  in  F.go  2 and  tnc 

• iw  Fiaure  5 shows  the  noise  reduction 
• Fin  i and  Fige  4,  respectively.  Figure  > sn; 

noisy  m Fig.  ) ana  riy.  , r 

d ui. j c filter  should  be  extended  by 
filter  used  in  synthetic  highs  system.  But  this 

. 3X3  window  so  as  to  pass  both  the  positive  and  the  negative  parts  of  boun- 
daries detected  by  Laplaclan  or  gradient  edge  detector,  otherwise  we  ose 

. • -nniflcantly  Reconstructed  pictures  are  given  ,n  Fig.  an 

resolution  s i gn i f i cant  I y^ 

• i maJian  filter  are  shown  in  Fig..  ° ana 
Fig.  7.  Outputs  of  the  two-dimensional  median 

, c r mpd i an  filtered  picture  we  lost  resolution  quite  a bit. 

Fig.  9-  In  5x5  median  filtered  P 

. effective  for  removing  impulse  errors  (all  dark 

Median  filters  seem  to  be  very  effective  . 

unpleasant  level.  Variable  width  Gaussian  filters  are  truncated  at  x 
standard  deviationo 

a eomif.ni- lal  lv  in  horizontal  ant  vertical 
A one-dimensional  filter  was  used  sequentially 

i ,k.  filter  in  the  vertical  direction  we  use  data 
directions  and  when  we  apply  the  f 
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already  averaged  in  the  horizontal  direction.  This  may  cause  more  reduction 
of  noise  but  resolution  Is  lost  also.  Actually  we  didn't  smooth  out  at  the 
top  interval  where  the  absolute  gradient  Is  highest.  (Filter  size  was  specified 
by  0 in  section  2.3.)  Comparing  with  the  variable  width  Gaussian  filter, 

the  performance  of  the  variable  width  median  filter  Is  very  poor.  In  the  noise 
cheating  technique  we  first  averaged  the  noisy  picture  by  2x2  window  and  re- 
placed the  grey  levels  at  4 picture  points  in  the  window  by  the  averaged  level. 

We  did  the  same  thing  with  3x3  window  and  el  Imlnated  Isolated  picture  blocks 
(Note:  3x3  window  will  have  the  same  grey  level  after  averaging)  by  simply 

replacing  the  center  grey  level  by  the  surrounding  grey  levels.  Since  the 
whole  block  (2x2. 3x3)  has  the  same  grey  level,  the  output  of  the  noise  cheating 
technique  has  lots  of  square  blocks.  The  performance  of  this  technique  seems 
the  worst.  The  two-dimensional  median  filter  is  most  economical  and  easiest 

to  apply  but  the  synthetic  highs  system  gives  the  best  result  although  It 
is  most  costly. 

Some  of  the  techniques  can  be  modified  for  better  performance.  For 

example,  we  may  use  the  noise  reduction  filter  used  in  the  synthetic  highs 

system  for  variable  width  Gaussian  or  median  filter  and  we  may  use  moving  3x3 

overlapping  window  replace  the  center  grey  level  by  the  averaged  level  rather 

than  assigning  the  same  grey  level  to  the  whole  block.  Finally,  we  emphasize 

that  the  performance  of  noise  reduction  techniques  depends  upon  the  type  of 
noises  involved. 
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Figure  k Low  pass  picture 
BW  - 0.116 


Figure  5 Noise  reduction  filter 
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Figure  6 Reconstructed  picture 
Laplacian  synthetic 
highs  used 


Figure  7 Reconstructed  picture 
Gradient  synthetic 
highs  used 


8 Reconstructed  picture 
3x3  median  filter  used 


Figure  9 Reconstructed  picture 
5x5  median  filter  used 
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Flgjure  III  Reconstructed  picture 
Variable  width)  Causslam 
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Figure  12  Reconstructed!  picture 
Variable  width  median 
filter  uni  form 
quantization  used 


Figure  13  Reconstructed  picture 
Variable  width  median 
filter  non-uniform 
quantization  used 
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Figure  14  Reconstructed  picture 

Noise  cheating  technique 
2x2  and  3x3  windows  used 
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TWO-DIHEMS IONAL  DIGITAL  RECURSIVE  FILTERS  AMD  THEIR 
APPLICATION  TO  IMAGE  PROCESSING 

Brian  O'Connor  and  T.S.  Huang 

Two  ft i mens i ona I digital  recursive  filters  have  the  potential  of  saving 
computer  time  and  storage  in  the  processing  of  large  two  dimensional  arrays  sucn 
as  images.  In  our  research  we  are  concerned  with  their  aesign  and  their  appli- 
cation to  images.  The  desired  processing  to  be  performed  or,  an  image  can  gen- 
erally be  described  in  the  frequency  domain.  In  two  dimensional  filter  design 
the  filter  coefficients  are  found  which  best  approximate  this  desireo  frequency 
response,  while  at  the  same  time  guaranteeing  stability. 

Several  possible  design  approaches  exist.  One  relies  on  nonlinear  optimi- 
zation techniques  to  minimize  the  t.  norm  of  the  difference  between  the  desired 
magnitude  or  group  delay  response  and  the  filter  response  [1,2],  These  tech- 
niques can  be  modified  to  guarantee  stability.  Hoy/ever,  numerical  problems 
arise  with  nonlinear  optimization  techniques,  namely  very  large  amounts  of 
computation,  sensitivity  to  starting  points,  and  the  possibility  of  converging 
to  a local  rather  than  global  minimum.  In  addition,  it  is  necessary  to  check 
stability  of  the  designed  filter  at  each  iteration  of  the  algorithm.  However, 
a recently  proposed  method  has  eliminated  need  for  the  stability  tests  [3J. 
Another  approach  is  based  on  an  algorithm  which  allows  the  designer  to  specify 
an  arbitrary  magnitude-squared  characteristic  which  is  approximated  optimally 
in  a weighted  Chebyshev  (minimax)  sense.  Here,  the  optimization  procedure  can 
be  formulated  as  a linear  programmi ng  problem  and  trie  filter  coefficients  can 
be  calculated  using  the  two-dimensional  discrete  Hilbert  transform  to  approxi- 
mately factor  the  two-dimensional  magnitude-square  frequency  response  [4], 

There  are  two  problems  in  using  this  algorithm.  The  first  is  that  fairly 
large  amounts  of  computer  time  are  needed  to  aesign  filters.  Secondly,  the 
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factorization  of  magnitude-square  function  to  obtain  a filter  implementation  con- 
tains an  infinite  number  of  terms  and  hence  truncation  is  necessary. 

Generally,  the  nonl inear  optimization  techniques  require  that  the  filter 
consist  of  first  and  second  order  sections  connected  in  cascade.  This  cascade 
form  has  many  advantages  over  the  direct  form  [t].  In  our  research  we  are  In- 
vestigating the  approximation  problem,  i.e.,  how  well  can  a cascade  of  first 
and  second  order  sections  approximate  an  arbitrary  frequency  response. 

Another  design  approach  is  to  find  filter  coefficients  which  approximate 
a desired  frequency  response  without  adding  the  stability  constraint;  then,  if 
the  filter  is  unstable,  stabilize  it  so  that  the  resulting  filter  has  approxi- 
mately the  same  frequency  response.  Several  stabilization  methods  exist.  The 
first  was  proposed  by  Shanks  [5l.  It  stabilizes  an  unstable  filter  by  finding 
its  double  planar  least  squares  inverse  (PLS I ) . Until  recently  the  PLSI  of  a 
filter  was  assumed  to  be  always  stable.  However,  a special  counter  example  has 
been  constructed  by  Genin  and  Kamp  where  they  find  a 2x2  PLSi  of  a 4x4  array 
[6].  It  is  still  an  open  question  whether  a PLSI  of  the  same  size  as  the  Input 
array  is  stable.  Jury  [7]  has  proved  this  for  the  special  cases  of  3x2  and  2x3. 
The  double  PLSi  of  an  array  will  produce  a filter  array  whose  magnitude  response 
approximates  that  of  the  unstable  filter.  However,  in  many  applications  the 
approximation  is  not  adequate,  in  our  work  we  have  found  an  example  where  the 
double  PLSI  calculated  by  Read  and  Treitei  is  unstable.  But  no  claims  on  find- 
ing a significant  counter  example  will  be  made  until  Read  and  Treitei's  [8] 
calculation  of  PLSI  can  be  checked. 

Another  approach  stabilizes  an  array  by  developing  a two-dimensional  dis- 
crete Hilbert  transform  to  calculate  the  analytic  phase  function  from  the  mag- 
nitude-squared frequency  response.  The  method  accomplishes  stabilization  with 
little  accompanying  distortion  of  its  amplitude  spectrum.  Several  developments 
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of  the  discrete  Hilbert  transform  exist*.  Read  and  Treitel  [8]  generalized  it 
to  two  dimensions  by  a direct  extension  of  the  one  dimensional  discrete  case. 
The  amplitude  spectrum  of  the  stabilized  filter  more  closely  approximates  the 
desired  magnitude  than  the  double  PLSI  filter.  However,  not  all  filters  ,->ro 
duced  by  this  method  are  stable.  An  example  of  this  was  given  in  their  paper 
and  we  have  found  many  more.  Another  method  relies  on  discretizing  the  con- 
tinuous two-dimensional  Hilbert  transform  to  obtain  a two-dimensional  discrete 
Hilbert  transform.  This  has  been  derived  by  Dudgeon  [9]  and  we  have  formulated 
and  programmed  a special  case.  Preliminary  results  seem  to  indicate  that  this 
method  does  produce  stable  filters,  but  the  approximated  amplitude  is  notable 
d i s tor ted . 

We  have  been  studying  a third  method  which  either  employs  the  cepstrum  or 
complex  cepstrum.  The  study  was  motivated  by  work  done  by  Ekstrom  and  Woods 
[10]  on  two  dimensional  spectral  factorization.  The  inverse  cepstrum  of  the 
first  quandrant  of  the  autocorrelated  unstable  filter's  cepstrum  is  windowed 
to  give  a stable  filter  whose  frequency  response  is  close  to  that  of  the  un- 
stable filter.  Preliminary  results  show  that  even  though  stability  is  not 
guaranteed  the  resulting  filters  are  usually  stable.  This  method  can  also  be 
used  to  test  stability  of  any  two-dimensional  recursive  filter.  The  original 
filter  is  stable  if  the  calculated  3rray  is  equal  to  the  original  array.  Be- 
cause the  FFT  is  used  the  correspondence  is  only  approximate  and  some  equality 
measurement  must  be  used  to  ascertain  stability.  Through  experimentation  we 
found  that  for  3x3  arrays  the  filter  is  stable  if  the  mean  square  difference 
between  input  and  output  arrays  is  less  than  .00000^3.  (16x16  FFT  were  used.) 

Another  possible  way  of  checking  stability  and  stabilizing  unstable 
filters  is  to  work  with  the  complex  [12]  cepstrum  of  the  filter  array.  By 
o roper ly  p-ocessing  the  cepstrum  we  can  obtain  arrays  which  are  nonzero 
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“rtai"  d'Stin9UiShed  r'S,°"5  «•  <V  V Pl«.  and  thus,  depending  on  the 

the  type  of  filter,  guaranteed  to  be  stability.  Many  problems  exist  in  one- 

dimenslnnal  cepstra,  analysis  and  they  are  increased  in  two  dimensions.  A two 

dimensional  complex  cepstrum  program  has  been  written  using  a modification  of 

a new  phase  unwrapping  algorithm  which  was  developed  by  Tribolet  [II].  Oet.il. 

about  the  theory  and  Implementation  of  two-dimensional  cepstrlal  analysis  will 
be  given  in  the  next  report. 

Future  worh  inciud.  a detailed  analysis  of  two-dimensional  cepstral 
techniques  applied  to  filter  design  and  stabilization. 

BIBLIOGRAPHY 

Oig^ta i*3 Recursive *FM ters"'"* ,EEE  Tra'"9"  for  T“°  Dimensional 

Feb.  1974.  r i Iters,  IEEE  TranSo  on  ASSP,  vol.  ASSP-22,  pp.  51-21 , 

ssi~-jsa,“  nj-jj: 

'51  ^i«"!{i^a!!d|S*{a.,,l,S!,,tyH,:f  Svnth'sls  0f  Dimensional 

Mters,  IEEE  Trans.  Audio  and  Electro.,  vol.  AU-20,  June  1972 

” mS^iS^**** 

[8]  R.  Read  and  S.  Treitel  "Th»  i i , _ , 

W s br* 

’9)  CaJr?d3ge0%aIrf97TSi0na'  ReC“rSiVe  Fi,teri"S."  Ph.O.  Thesis,  MIT, 


• v,7.VW*^rt-  f'CTS1  rX1 


69 


[10]  M„  Ekstrom  and  J„  Woods,  "Two-Dimensional  Spectral  Factorization  with 

App  I i ca  t i ons  In  Recursive  Digital  Filtering,"  IEEE  Trans.  ASSP,  vol.-2l4, 
pp.  1 15-127,  Apri I 1976- 

fl!]  J.  Iribolot,  "A  New  Phase  Unwrapping  Algorithm,"  submitted  to  IEEE  Trans, 
on  Acoustics,  Speech  and  Signal  Processing,  March  1976. 


[12]  A.  Oppenheim  and  R.  Schafer,  Digital  Signal  Processing,  Prentice-Hall, 

Inc.  1975. 


70 


COMPARISON  OF  THE  PROJECTION  METHOD  WITH  SINGULAR 
VALUE  DECOMPOSITION 

S.P.  Berger  and  T.S.  Huang 

The  purpose  of  the  work  has  been  to  evaluate  the  relative  merits  of  the 
projection  method  of  image  restoration  as  compared  with  the  singular  value  de- 
composition (SVD)  approach. 

The  projection  algorithm  is  an  iterative  method  of  solving  a set  of  linear 
equations,  where  the  equations  are  the  discrete  representation  of  the  degrada- 
tion process.  The  algorithm  can  incorporate  a priori  information.  It  has  been 
shown  to  yield  effective  results  for  various  types  of  degradation,  including 
space-variant  distortions. 

The  SVD  approach  has  been  used  successfully  in  image  restoration.  It  in- 
volves the  treatment  of  the  degradation  in  matrix  form.  By  application  of  a 
pseudo-inverse  matrix,  the  restoration  is  hopefully  achieved. 

Both  approaches  are  limited  by  the  effects  of  noise.  The  required  number 
of  iterations  in  the  projection  method  and  the  number  of  terms  utilized  with 
SVD,  are  determined  subjectively.  The  effects  of  noise  increase  with  the 
number  of  terms  and  the  number  of  iterations,  and  tend  to  overshadow  the 
restoration  process. 

In  the  actual  implementation  of  the  SVD,  difficulties  with  this  approach 
have  arisen.  The  amount  of  computer  time  required  for  the  calculation  of 
eigenvalues  and  eigenvectors  can  be  prohibitive  for  large  degradation  matrices. 
Also,  perhaps  to  the  size  of  the  matrix,  the  actual  implementation  has  yielded 
faulty  results.  The  success  of  the  method  depends  heavily  on  the  type  of 
degradation  that  is  effected  on  the  border  of  the  image.  Minor  changes  in  the 
form  of  the  degradation  matrix  seem  to  create  major  problems  in  the  operation 
of  the  computer  implementation.  The  next  report  will  contain  at  least  a 
partial  resolution  of  these  difficulties. 
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FOURIER  DESCRIPTORS 
T.  Wa 1 1 ace  and  P.A.  Wlntz 

The  Fourier  descriptor  (FO)  is  one  method  of  describing  the  shape  of  a 

i ni--~  the  contour  can  be  traced, 
planar  figure.  Given  a figure  in  the  comp  ex  p . 

yielding  a c»p.ex  function  of  t,«.  .»  - -our  is  traced  repeated 

periodic  function  which  r.suits  can  he  expressed  in  a Fourier  series. 

[0  defines  the  FO  of  a contour  as  the  coefficients  of  this  Fourier  series 

TP  implement  this  method  of  shape  description,  it  is  necessary  to  sample 

* a seguenc.  gives  us  the  values  of  the  Fourier  series  coefficients  of  the 
seguence,  assuming  It  to  he  periodic,  using  an  FFT  algorithm  s.tis  es  • 
definition  above.  The  cputatlona.  advantages  of  the  FFT  are  w.  no  . 

The  goal  of  this  worh  Is  to  classify  the  shapes  of  objects  using  the.r 

Fourier  descriptors.  The  operations  of  rotation,  scaling,  and 

starting  point  are  easily  implemented  In  the  freguency  domain  by  s,mp  « 

me,,  on  the  freguency  detain  coefficients.  Wh„.  shapes  may  be  ^ 

the  space  domain,  the  procedures  reguired  to  adjust  their  sit.  an  o 

ii  ifaraflue  type  of  algorithm  Is 
T ii..  expensive.  Normally  an  iterative  cyp 

are  computationally  very  expensive. 

tM‘  ::,o.l  of  class, f.catio Fourier  descriptors  is  to  develop  an 

eigorithm  which  will  normalize  the  size  and  orientation  of  a shape  before  any 

comparisons  to  test  shapes  are  made.  If  this  can  be  accomplished,  the  Cass  - 

•mole  clustering  problem  with  no  iterative  searches 
fication  process  becomes  a simple  clustering  p 

to  contend  with. 

, , the  FD  normalization  problem  was  dis- 

ln  our  last  quarterly  report,  th 

j i f wss  shown  *th3t  i t tine 
cussed,  and  the  results  of  an  experiment  presente  . 
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contours  under  investigation  were  bilaterally  symmetric,  much  less  Information 
was  contained  in  the  phases  of  the  FO  than  In  the  magnitudes.  The  validity  of 
classifying  contours  using  the  magnitudes  of  their  FD's  was  demonstrated  by 
classifying  18  airplane  contours  using  this  method.  The  Euclidean  distance 
measure  was  used,  and  the  distances  between  FDs  proved  to  correlate  well  with 
the  actual  differences  between  the  contours. 

The  goal  of  this  study  of  Fourier  descriptors  Is  to  eventually  apply  this 
method  to  contours  traced  using  actual  photographic  data.  To  simulate  this 
situation,  an  experiment  was  conducted  using  20  aircraft  contours  digitized  to 
two  different  resolutions,  128x128,  and  64x64.  The  high  resolution  versions 
were  taken  to  be  accurate  representations,  and  the  lower  resolution  versions 
were  assumed  to  be  corrupted  by  noise  and  quantization  error.  Examination  of 
representative  contours  (Figs.  I —6)  show  that  the  128x128  contours  are  quite 
good  representations,  while  the  64x64  show  significant  distortion  of  the 
smaller  important  features.  This  experiment  was  performed  In  order  to  test 
various  distance  measures,  as  well  as  to  test  suitability  of  the  algorithm  for 
use  with  actual  photographic  data. 

While  the  experiment  described  in  [1]  was  useful  In  demonstrating  the 
general  validity  of  classifying  contours  using  distances  between  normalized 
Fourier  descriptors  (NFDS),  a comparison  of  various  distance  measures  was 
difficult  due  to  lack  of  a definitive  measure  of  similarity  among  the  contours 
themselves.  This  obstacle  was  overcome  by  using  the  different  resolution 
versions  of  the  same  planes. 

The  mean  square  criterion  used  previously  was  compared  to  an  absolute 
value  criterion.  The  results  using  the  absolute  value  criterion  were  slightly 
better,  as  every  64x64  contour  was  correctly  identified,  whereas  using  the 
mean  square  criterion,  only  19  out  of  20  were  correctly  identified. 
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A simple  application  of  Parseval's  theorem  shows  that  the  mean  square 
criterion  in  the  frequency  domain  corresponds  to  a sum  of  the  point  by  point 
mean  square  error  in  the  time  domain.  There  are  some  difficulties  in  mathe- 
matically relating  the  frequency  domain  absolute  value  criterion  to  the  con- 
tours in  the  time  domain.  However,  it  is  easy  to  compare  the  two  criteria  on 
a qual i tat i ve  bas i s. 

It  is  obvious  that  the  absolute  value  measure  will  be  more  tolerant  of 
one  or  two  large  coefficient  differences  than  the  mean  square  measure.  Using 
either  method,  the  largest  coefficients  will  account  for  most  of  the  distance 
measured.  Accordingly,  we  should  consider  the  possibility  of  large  variations 
between  coefficients  of  large  expected  value,  which  turn  out  to  be  the  lowest 
frequency  ones. 

The  second  largest  coefficient  for  each  airplane  TD  is  A(-3).  The  effects 
of  varying  this  coefficient  are  to  change  the  width  of  both  the  wings  and  the 
body  of  the  plane.  Smaller  detail  such  as  engine  shape  is  virtually  unaffect- 
ed. Another  coefficient  of  large  expected  value  is  A(-l).  As  discussed  in 
[1],  this  coefficient  describes  the  length  to  wingspan  ratio  of  an  airplane 
contour.  More  generally,  varying  A (— l ) tends  to  elongate  any  contour.  Again, 
the  smaller  detail  is  not  greatly  changed. 

In  summary,  the  absolute  value  criterion  should  tolerate  more  differences 
in  gross  structure  of  the  contour,  such  as  elongatedness  or  thickness  to 
length  ratio,  while  emphasizing  variations  in  smaller  detail.  It  seems  likely 
that  the  absolute  value  measure  might  correspond  more  clearly  to  the  differ- 
ences which  human  observers  find  important  than  does  the  mean  square  criterion. 

The  classifications  using  the  two  methods  were  too  similar  to  offer  a 
definitive  comparison  of  the  two  classification  methods.  However,  the  success 
of  the  experiment  shows  that  we  are  ready  to  apply  the  FD  algorithm  to  contours 
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extracted  from  actual  photographic  data.  We  plan  to  interface  the  BLOB  con- 
tour tracing  algorithm  to  the  FD  program  as  the  next  step  in  this  research. 
The  magnitude  vs.  phase  information  question  will  also  be  examined  in  the 
context  of  developing  a normalization  procedure  which  preserves  all  of  the 
information  contained  in  the  contour. 
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filtering  to  remove  cloud  cover  in  satellite  imagery 

0.  R.  Mitchell  and  P.  L.  Chen 

I.  INTRODUCTION 

This  recent  work  is  a continuation  of  the  previous  report  "Satellite 
Imagery  Noise  Removal"  [1],  We  are  using  3-dlmenslonal  homomorphic  filtering 
techniques  to  remove  cloud  cover  In  LANDSAT  data. 

In  the  previous  filtering  results  the  noise  power  spectrum  was  estimated 
by  classifying  each  region  of  the  noisy  picture  according  to  the  level  of  noise 
present  using  the  mul tispectral  data  analysis  software  system  developed  by 
Lars.  These  filtering  results  were  not  satisfactory  due  to  the  spectral 
indlstinctabi 1 ity  of  clouds  and  concretes. 

it  may  also  prove  Impractical  to  use  a mul tispectral  classification  pro- 
gram to  find  the  noise  statistics.  Instead  it  may  be  possible  to  use  a 

generalized  cloud  power  spectrum  derived  by  averaging  many  sample  spectrums 
together. 

I • . INDIRECT  ESTIMATION  OF  NOISE  STATISTICS 

• he  easiert  way  to  estimate  the  general  power  spectrum  of  cloud  Is  done 
by  using  data  taken  over  water  where  the  reflection  Is  almost  constant,  in 
this  case,  the  transformed  scanner  Image  (L  Is  sun  Illumination,  t Is  cloud 
transmission,  r Is  ground  reflection,  and  s Is  the  received  scanner  Image) 

Log  [L  - s(x,y) J - Log  L + Log  t(x,y)  + Log  [L  - ar(x,y)] 
is  reduced  to 

Log  [L  - s(x,y)]  « Log  [L  - ar(x,y)]  + K 

where  K is  a constant,  and  the  power  spectrum  obtained  Is  that  of  noise  except 
for  the  d.c.  (0,0)  frequency  point. 
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This  power  spectrum  should  be  circularly  symmetric  since  clouds  have  no 
preferred  orientation  and  should  consist  mainly  of  low  spatial  frequency  com- 
ponents since  clouds  are  relatively  large  and  smooth  functions  compared  to 
ground  reflectance. 

III.  THREE  DIMENSIONAL  FILTERING 

The  raai  potential  In  the  cloud  filtering  process  is  Incorporating  a 
third  dimension,  the  spectral  channels,  forming  a three  dimensional  reflection 
r(x,y,z)  and  cloud  transmission  t(x,y,z).  The  generalized  linear  filter  thus 
employed  Is  three  dimensional  H(p,v,p)  using  three  frequencies  (two  spatial 
and  one  spectral).  Although  there  are  only  four  points  In  the  spectral 
dimension  for  LANDSAT  data,  the  method  has  good  premise,  because  most  clouds 
follow  a fixed  response  in  the  spectral  dimension:  cloud  transmission 

Increases  with  wavelength  In  a predictable  fashion.  When  this  Information  Is 
Incorporated  Into  the  filter  (by  means  of  the  3"D  power  spectrum)  image 
variations  which  have  the  cloud  spectral  response  are  filtered  out  and  image 
variations  which  do  not  follow  ' s expected  response  of  clouds  in  the  spectral 
dimension  are  left  In.  The  three  dimensional  filter,  therefore,  tends  to 
reject  all  variations  that  are  low  frequency  In  the  spatial  dimension  and 
follow  the  cloud  spectral  response  In  the  third  dimension. 

Figures  1 to  4 are  the  3~D  power  spectrum  obtained  from  averaging  the 
spectrum  from  three  separate  64x64  regions  of  clouds  over  water.  The  ordinate 
is  a log  scale.  Figure  1 represents  the  spectral  d.c.  slice  (sum  of  all 
four  spectral  points).  Fig.  2 represents  the  "one  cycle  per  spectral  band" 
slice,  etc.  One  Interesting  observation  we  have  made  is  that  most  of  the 
cloud  information  is  contained  In  slice  I only  and  that  information  in  slices 
2,  3.  and  4 (other  than  the  d.c.  point;  in  each)  represent  ground  reflectance 
Information  and  should  be  left  in  the  filtered  output. 
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Figures  5 to  8 are  the  3-0  filter  functions  obtained  for  one  particular 
ikxbk  section  of  the  Image  taken  over  land.  The  general  cloud  spectrum  was 
normalized  so  that  it  fell  below  the  spectrum  of  "signal  plus  noise"  and 
was  then  subtracted  from  It  for  the  signal  spectrum  estimate.  Note  the  filter 
tends  to  attenuate  low  frequencies  In  slice  one  but  leaves  them  In  slices  2 
through  k. 

Figure  9 shows  the  composite  filtered  for  21  6i*x6i*  blocks  of  LANDSAT  data 
scanned  over  the  Chicago  area  on  a somewhat  cloudy  day. 

IV.  CONCLUSIONS 

Three  dimensional  filtering  of  multispectral  data  to  remove  light  cloud 
cover  is  a distinct  possibility. 

In  Fig.  9 the  filtered  picture  shows  more  detail  than  the  original  cloudy 
one.  Lake  boundaries  and  highways  have  clearer  appearances  In  the  filtered 
pictures. 

The  model  of  the  cloud  distortion  needs  to  be  refined  based  on  the 
results  of  filtering  using  the  simple  model  presented  In  the  previous  report. 
It  may  be  necessary  to  consider  convolutional  effects  of  cloud  cover  as  well 
as  multiplicative  effects.  The  change  In  multispectral  classification 
accuracy  after  filtering  may  be  a suitable  measure  of  the  performance  of  such 
homomorphic  filter. 
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Figure  2 3-0  cloud 
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Figure  3 3-D  cloud  power  spectrum 
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Figure  5 3”D  filter  function  - slice 


Figure  6 3-0  filter  function  - slice 


Figure  7 3-D  filter  function  - slice 
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Figure  8 3-D  filter  function  - slice 
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Figure  9(a)  (Left)  3-D  filtered  output  (channel  3) 

(Right)  Original  Landsat  data  (channel  3) 


Figure  9(b)  (Left)  3-D  filtered  output  (channel  **) 

(Right)  Original  Landsat  data  (channel  **) 
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3 

7 

1 


Manufacturer 
Beehive  Elect. 
Tex.  Inst. 
Dlgl-Oata 


DEC 

DEC 

Fabrl tek 

Versa tek 
Comtal 

Data  Printer 
True-Data 
Tektronix 
DEC 


FACILITIES 


Description 
"Super-Bee"  Terminals 


"Silent  700"  Terminals 

Industry  standard  magnetic  tape  system; 

2,  9-track  and  1,  7-track  drives;. 

NRZI  and  phase-encoded  formatters/control le 


Dual -drive  DEC tape  unit 

RP03  disk  drive  (40  million  characters) 

96K-word  auxiliary  memory  system 
(64K  bought  by  ARPA,  32K  by  NASA) 

Electrostatic  matrix  printer 


Color  picture  display 

132  column,  600  L.P.M.  line  printer 


Punched  card  reader 

Model  4010,  graphics  display 

PDP-11/45  computer  system;  system  Includes: 

32K  memory  . _ * 

FPP-11  floating  point  processor  (NSF  money) 

H960  extension  mounting  cabinet  . 

3 - small  peripheral  mountings  blocks  (DD  11) 

1 UNI  BUS  repeater/expander 

DH1 I , 16-1 tne  terminal  multiplexor 

KWll-p  programmed  clock 

"ANTS"  - type  PDP-11 /IMP  Interface 


Note:  Our  PDP-ll/45  Is  currently  operating  under  the  UNIX  system 
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